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During this period investigations were carried out in three areas: 


1. Binding energies and high energy binding sites were calculated for 

carbon atoms deposited on a (2x1) dimerized Si(100) surface. The 
location of the energetically most favorable binding site was found 
to be near the top position of a second layer Si atom with a bind- 
ing energy of —3.36 eV. Also calculated were excess energies for 
the substitution of C (replacing Si atoms) in the exposed surface 
region. For the substitution of a carbon atom in surface layers 1 
through 4 calculations produced progressively increasing excess 
energies from the exposed surface to the interior. This indicates 
that the diffusion of a substituting C atom from the top layer to 
the interior via a site exchange mechanism is energetically unfa- 
vorable. Computational details along with a short discussion of 
the results are presented in Appendix I. 

2. Structure- and energy-related properties of defects formed on re- 

laxed diamond surfaces were investigated. Simulation calculations 
were carried out for vacancies created on (ill) and on (2x1) dimer- 
ized (100) index planes of diamond. Two different model functions 
(based on two- and three-body interactions) developed recently for 
carbon were employed in the calculations. Both functions produced 
comparable results. Top layer vacancies are more likely to form on 
the (100) surface than on the (111) surface. Calculations indicate 
that the diffusion of a single vacancy from the top surface layer 
to the second layer is not energetically favored. For details see 
Appendix II. 
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3. In this part, recently developed Tersoff-type model functions for Si 
were employed to investigate the energetics of the following pro- 
cesses: (1) neutral monovacancy and divacancy formation and migra- 
tion energies; (2) neutral bond-centered, site sharing, tetrahe- 
dral, and hexagonal self -interstitial formation and migration en- 
ergies; and (3) the variation of these energies with distance from 
a bulk site of Frenkel defect formation and from [100] and [111] 
surfaces. All the calculations were performed employing an energy 
minimization procedure for the low temperature limit. In addition 
to the energy values, structural changes around defective regions 
were also calculated for different vacancies and interstitials. 
Details of this investigation are presented in Appendix III. 
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APPENDIX I: 


Carbon Atoms Deposited on the (2x1) Reconstructed Si(100) Surface 
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Introduction 





An atomic level understanding of the structure of a surface defect 
and its energetics is highly desired by many surface scientists today 
[1] . In diamond film synthesis, in particular, it has been shown that 
defects in the surface region of a substrate play an important role in 
early stages of the nucleation and growth process [2,3] . Binding ener- 
gies for carbon atoms on Si substrate surfaces along with the binding 
sites geometries are very important features which often determine the 
quality and characteristics of the film. 

In this investigation simulation calculations were carried out to 
analyze energy- and structure-related properties for carbon atoms de- 
posited on a (2xi) dimerized Si(100) surface. In addition to binding 
energies, substitutional energies were also calculated for C atoms re- 
placing Si atoms in the surface region. All the calculations in this 
study were carried out for the low temperature limit using a static min- 
imization procedure. Interatomic energies as functions of atomic po- 
sitions were estimated employing a potential energy function developed 
recently by Tersoff [4] . This function has been designed specifically 
for SiC and it has been used successfully in calculating energies and 
structures for carbon defects in bulk Si [5] . At the same time, it has 
been shown that for pure Si systems also, the Tersoff function produces 
acceptable results for various bulk and surface properties [6] . Fur- 
thermore, for pure carbon systems this function is able to reproduce 
correctly many bulk properties of diamond and of the graphitic plane, 
and recently it has been used to estimate surface properties for low in- 
dex planes of diamond [7] . 
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Calculations 
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The Si (100) surface was generated, first, as an abrupt termination 
of bulk silicon in the diamond cubic structure with a lattice constant 
corresponding to the equilibrium volume of Si calculated by the Tersoff 
function. Then, the system was equilibrated producing (2x1) dimerized 
reconstruction patterns at the exposed (100) surface. In the equilibra- 
tion process, every atom in the system was permitted to relax by mini- 
mizing the total energy with respect to the atomic coordinates. Simu- 
lation calculations were carried out considering a computational cell 
which was made of a slab of 16 atomic layers each containing 16 Si atoms. 
This cell size was found to be adequate for the present investigation. 
Calculations for selected cases using larger cell sizes produced virtu- 
ally identical results. Throughout this study, calculations were per- 
formed considering periodic boundary conditions imposed on the system 
in two directions (parallel to the exposed surface) in order to provide 
continuity. 

The binding energy, <f>b, per carbon adatom was calculated as: 

<f>b = -[ 4 m) - E°] 

m 

where m is the number of adatoms, E° denotes the total equilibrated 
energy of the system of N particles with clean exposed surfaces, and 
E is the total relaxed energy of the same system with m adatoms de- 
posited on the surface. In calculating binding energies, carbon atoms 
were first positioned on a fully relaxed surface and then the system was 
reequilibrated to obtain E ^ 1 ^ . In order to find adsorption sites with 
lowest binding energies, calculations were repeated with different ini- 
tial positions for carbon atom. On a single surface cell 25 different 
initial configurations were employed and those producing lowest ener- 
gies were reported here. 
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In addition to binding energies, we also calculated energies of sub- 
stitution for C replacing Si atoms in the surface region. Substitutions 
were made in different surface layers and the system was fully reequil- 
brated after the substitution. In a manner analogous to the binding en- 
ergy, the substitutional energy also was considered simply as an excess 
energy in this investigation and it was calculated as: 


= i [£<">- E”) 

n 

where (f) x represents substitutional excess energy, n is the number of 
substituting C atoms, E° denotes the total equilibrated energy of the 
Si system bearing the exposed surface (as defined above) and E s is the 
total relaxed energy of the same system with n substitutional C atoms. 

In the process of substitution the number of atoms in the system remains 
unchanged. Furthermore, according to the present definition <f> x repre- 
sents an excess energy, therefore, it does not include formation ener- 
gies for pure C and Si systems. 

Results and Discusion 

For a single carbon atom the most favorable binding site on a defect- 
free (2x1) dimerized Si (100) surface was found to be ‘almost’ on top of 
a second layer Si atom between two dimer rows. Schematic top and side 
views are shown in Figure 1. Due to strong Si— C interactions Si atoms 
located in the vicinty of C were displaced from their original sites. 
Upon relaxation all neighboring Si atoms had migrated toward the C atom. 
In this case, there are three neighboring Si atoms in close proximity to 
the deposited C atom. The nearest one is located in the second layer at a 
distance of 1.72 A. The other two neighbors are in the top surface layer 
and belong to two different Si dimers. They are equidistant from the C 


6 



adatom and the Si — C bonds were calculated as 1.82 A forming an apex an- 
gle of 133 degrees. The binding energy for the carbon atom at this site 
was calculated as —6.193 eV, 

The other energetically favorable binding site for a single C atom 
is located above the mid-point of a Si dimer. This binding site is Ex- 
actly’ above a fourth layer Si atom and it is depicted in Figure 2. In 
this case, the C adatom has two close neighbors with equal Si— C sepa- 
rations of 1.80 A. These bonds form an apex angle of 124.5 degrees at 
the carbon atom position. Under the force field of the carbon atom, 
the neighboring Si atoms forming the dimer were displaced toward their 
original (lxl) positions. Here, in its fully relaxed configuration the 
Si— Si separation was calculated as 3.18 A, while for the defect-free 
(2x1) reconstructed Si(100) surface, the Si— Si distance in a dimer is 
2.37 A. The binding energy for a C atom at this site is somewhat higher 
than the previous case and it was estimated as —5.925 eV. 

For two adatoms the energetically most favorable sites are shown in 
Figure 3. These sites are symmetrically situated along the dimer row 
and they are basically identical with the high energy position indi- 
cated above for the single C atom. The four Si— C distances between the 
C adatoms and their proximate top layer neighbors in the adjacent Si 
dimers were calculated as 1.79 A. Distances between the carbon atoms and 
their nearest neighbors in the second layers were found to be equal to 
1.70 A. These distances are sligthly shorter than the single carbon atom 
case mentioned above. The total binding energy per C adatom is —6.30 eV 
which is somewhat lower than the energy value obtained for a single C 
atom calculated above. Because of the local symmetry in this case, an 
additional strain energy was released upon relaxation. The distance be- 
tween two carbon adatoms is 3.70 A and obviously at this separation no 
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contribution is expected coming from the C — C interactions to the total 
energy . 

Next, excess energies for substitution were calculated for a C atom 
replacing Si within the exposed surface region. Calculations were per- 
formed for the (2x1) dimerized (100) surface. Excess energies were es- 
timated for C atom substitutions in the surface layers 1 through 4. Cal- 
culated energy values along with some structural parameters are pre- 
sented in Table 1. 

An equilibrated top and side views are depicted schematically in 
Figure 4 for a single C atom replacing a top layer Si atom. In this case, 
one of the Si atoms in a dimer is replaced by the carbon atom. In the 
fully equilibrated configuration the site for the substitutional C atom 
is about 0.45 A lower than the top surface layer and the Si— C bond is po- 
sitioned asymmetrically. See Figure 4. In this dimerized top layer, the 
substitutional C atom has three neighbors with approximately equal Si— C 
bond lengths. While one of the Si atoms is in the first layer, the other 
two are positioned in the second layer. In this case, due to the relax- 
ation in the interlayer spacings, the avarage Si— C— Si angles were found 
to be around 119 degrees, somewhat larger than the ideal tetrahedral an- 
gle. 


Each C atom substituting for a Si atom in layers 2 through 4, has 
four immediate neighbors located in the contiguous layers. Two of them 
axe in the next upper layer and the other two are located in the layer 
below. For a C atom substitution in the second layer the excess energy 
was found to be somewhat higher than in the previous case. The average 
Si— C distances to immediate neighbors are slightly larger. (See Table 
1) . The Si— C— Si angles which can be regarded as a measure of the tetra- 
hedral symmetry in this case vary within 104 - 113 degrees averaging to a 
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value of 109 degrees. 


On a clean (2x1) dimerized (100) surface all the atomic sites in 
the top layer are identical. Similarly, in the second layer also all 
the sites are identical and therefore, they would produce the same ex- 
cess energy upon a C atom substitution. For the third and fourth layers, 
however, two different substitutional sites are possible. Figure 5 de- 
picts a schematic top view of a fully relaxed Si (100) surface. Atomic 
sites up to fourth layer from the top are shown in this figure. The two 
different substitutional sites for carbon in the third and fourth lay- 
ers were indicated in Figure 5 by letters <z, b and c, d, respectively. 
Substitutions in sites o and c which are located under the dimer row, 
produced lower excess energies than the substitutions in sites b and d. 
The substitution in site a produced a value comparable with the second 
layer substitutional excess energy obtained above. The value for site 6, 
however, is about 0.85 eV higher. These excess energies give an average 
(j> x value of —2.04 eV for the third layer substitution. A similar trend 
was obtained for the fourth layer as well. The site c, which is located 
beneath the mid point of a Si dimer, produced a <f> x value which is about 
0.76 eV lower then the site d. In this case, an average value of —1.87 
eV was obtained for the substitution of C in layer 4. Average Si— C bond 
lengths for both the third and fourth layers, however, show only small 
variations. Similarly, changes in the Si— C— Si angles were also found 
to be very small and calculations for the third and fourth layer cases 
produced average tetrahedral angles around 109 degrees. 

Conclusions 

Calculations carried out in this investigation indicate that the 
deposition of C atoms on a (2x1) dimerized Si (100) surface is energet- 



ically quit© favorable. For a single C adatom the lowest energy bind- 
ing site is located near the top position of a second layer Si atom be- 
tween two Si dimer rows. In this case the local symmetry at the binding 
site was somewhat distorted and the adjacent Si dimers were tilted to- 
ward the C atom. Substitutions of carbon atoms in the top surface layers 
were also found to be energetically favorable. Due to the strong Si— C 
interactions after the substitution the local symmetry of the lattice 
is distorted and in all cases, the neighboring Si atoms were found to be 
displaced toward the substituting C atom. 

Replacements of a Si atom by C in surface layers 1 through 4 produced 
progessively increasing average excess energies from the exposed sur- 
face toward interior layers. Accordingly, the diffusion of a substitut- 
ing C atom from the top layer to the interior via a site exchange mecha- 
nism is energetically unfavorable. Results obtained in this investiga- 
tion, however, are based strictly on energetics and no entropic aspects 
were taken into consideration. Therefore, it is recommended that ex- 
treme care be exercised when comparing these results with experimental 
findings . 
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Table 1. Calculated values for the replacements of Si atoms by C in the 
top four layers of the Si (100) surface. The excess energy for C substi- 
tution, the average Si— C distance and the number of near neighbors are 
denoted by <f> x , dsi-c and respectively. For the third and fourth lay- 
ers different substitutional sites indicated by letters a through d are 
shown schematically in Figure 5. (For first and second layers all sites 
are identical) . 


Layer 

Number 

Site 

4>x 

(eV) 

ds i-c 
(A) 

Ni 

1 

0 

-3.36 

1.85 

3 

2 

0 

-2.45 

1.91 

4 

3 

(a) 

-2.47 

1.93 

4 

3 

(b) 

-1.61 

1.95 

4 

4 

(c) 

-2.25 

1.94 

4 

4 

(d) 

-1.49 

1.97 

4 
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Figure Captions 
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Figure 1. Schematic (a) top and (6) side views showing the energetically 
most favorable site for a carbon adatom deposited on a (2 x 1) recon- 
structed surface of Si(100). Si atoms in dimers axe connected with solid 
lines and dotted lines indicate bonds between the top layer atoms and 
atoms located in the second layer. Large, medium and small open circles 
indicate sites for Si atoms located in the first, second and third lay- 
ers, respectively. The deposited carbon atom is represented by a large 
shaded circle. 

Figure 2. Schematic (a) top and (6) side views for a carbon atom de- 
posited on the Si (100) -(2 x 1) surface. The carbon adatom, in this case, 
is located above the mid-point of a Si dimer. (Symbols are explained in 
the caption for Figure 1 . ) 

Figure 3. Two carbon atoms deposited on a Si (100) -(2 x 1) surface. Top 
and side views are depicted in a and b> respectively. (Symbols are ex- 
plained in the caption for Figure 1.) 

Figure 4. In the top layer of a Si (100) -(2 x 1) surface a single C atom 
is shown substituting for a Si atom. Top and side views are depicted in a 
and by respectively. (For symbols see the caption of Figure 1). 

Figure 5. A schematic top view for a fully relaxed (2 x 1) dimerized 
Si (100) surface. Large, medium and small open circles indicate Si atoms 
located in the first, second and third layers, respectively. Small 
solid circles represent Si atoms in the fourth layer. Letters a, b and c, 
d indicate different sites for Si atoms located in the third and fourth 
layers, respectively. 
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Calculations for Defects Formed on Diamond Surfaces 
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Introduction 


Surface defects play a significant role in research areas related 
to nucleation and growth phenomena. An understanding of defects at mi- 
croscopic dimensions is highly desired by surface scientists today. In 
general, it is believed that the local structure of a surface defect is 
an important factor controlling nucleation and the early stages of the 
growth process [1,2] . Presently, there is not enough information avail- 
able on the energetics and the local structure of defects formed on dia- 
mond surfaces. Computer simulation techniques based on atomic consider- 
ations provide a very useful approach to the study of various processes 
taking place on surfaces . 

In this work, a simulation study is conducted to investigate 
structure- and energy-related properties of defects formed on the re- 
laxed (111) surface and on the (2x1) dimerized (100) surface of dia- 
mond. Simulation calculations were performed for the low temperature 
limit (i.e. , T=0 K) using a static minimization procedure. To calcu- 
late configurational energies model potential functions developed re- 
cently for carbon species were used throughout this study. For a better 
analysis, two different model functions, based on two and three-body in- 
teractions, were considered here. The first potential is the Tersoff 
function [3] . It has been shown by Tersoff that this function is able 
to reproduce correctly various bulk properties of diamond and of the 
graphitic plane. The second potential used in this study is the Brenner 
function [4] which is analytically similar to the Tersoff function, but 
parametrized differently. This function also, has been shown to produce 
acceptable results for properties of bulk diamond and the basal plane of 
graphite, as well as for some properties of small carbon clusters [4,5] . 
While the Tersoff function reproduces correctly the lattice constant 



of diamond, the Brenner potential at its minimum energy configuration 
gives a lattice constant value about 3*/, shorter. However, structural 
properties for small clusters of carbon are better represented by the 
Brenner function [5] . These functions have also been used recently to 
calculate defect-free surface properties for low index planes of dia- 
mond [6] . 


Method of Calculation 

In this study structure- and energy-related properties of small va- 
cancies created on the (111) and (100) surfaces of diamond were investi- 
gated. Calculations were carried out considering a computational cell 
containing more than 500 carbon atoms in a diamond cubic structure. The 
system was positioned in two different ways to provide the desired (111) 
and (100) index planes with layers parallel to the exposed surface, each 
containing 32 carbon atoms. Periodic boundary conditions were imposed 
on the system in two directions (parallel to the exposed surface) in or- 
der to provide continuity. These systems then were fully equilibrated 
to obtain relaxed surfaces. For the (100) surface the equilibration 
was carried out to produce a (2x1) dimerized reconstruction pattern. 
Defects were created on these relaxed surfaces and then the system was 
reequilibrated to obtain relaxed defect energies and geometries. 

In all cases, defect energies were calculated as: 

E v = - $ 0 + (Ni - N 2 ) x e coh 

where, E v is the surface vacancy formation energy, denotes the total 
relaxed energy of the defect-free system with an exposed surface. is 
the total relaxed energy of the same system with surface vacancies. e co h 
represents the cohesive energy of bulk diamond calculated by the poten- 
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tial function. Ni and N 2 denote respectively the number of carbon atoms 
in the system before and after the formation of the defect. To analyze 
the extent of the relaxation taking place during the reequilibration 
process we also calculated the unrelaxed energy, E^ , for defects. In 
this case, E ^ is calculated using the unrelaxed value of before the 
reequilibration. Accordingly the relaxation energy for the vacancy is 
defined as the difference between E Y and E v . 


Results 

For the (2x1) dimerized (100) surface, vacancy energies calculated 
using Tersoff said Brenner functions sore given in Table 1. A schematic 
top view of the dimerized (100) surface and vsirious defect structures 
are shown in Figure 1. For a single vacancy created at the top layer, 
both the Tersoff and Brenner functions produce relatively low forma- 
tion energies. The value calculated by the Tersoff function is some- 
what lower than the value obtained by the Brenner potential. In both 
cases, however, the relaxation lowered vacancy energies considerably. 
For the dimerized (100) surface, upon creating a single vacancy at the 
top layer, one carbon atom is left unpaired. During the reequilibra- 
tion process it was found that this unpaired C atom moves to its original 
(1 xl) position, which is right above the mid-point between two second 
layer C atoms (see Figure lb). 

Both potential functions produced much laxger energy values for a 
single vacancy created in the second layer. In this case, the Brenner 
function predicts somewhat lower energy than the Tersoff function and 
relaxation energies associated with the vacancies are quite large. Af- 
ter the reequilibration, C atoms adjacent to a second-layer vacancy were 
found to be displaced radially outward. The two C atoms closest to the 
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vacant sit© are located in the top layer. During the reequilibration 
process these two carbon atoms migrated outwards (in the direction away 
from the vacancy center) about 0.25 and 0.32 A for the Tersoff and Bren- 
ner functions, respectively. The top view of this relaxed configuration 
is shown schematically in Figure lc. 

A dimer vacancy was generated in the top layer by removing both car- 
bon atoms of a dimer. This is shown in Figure Id. Calculations for a 
dimer vacancy produced very low energy values. While the Tersoff func- 
tion produced a dimer vacancy energy comparable to the energy of a sin- 
gle vacancy, the Brenner function predicts even a lower energy value for 
a dimer vacancy. (See Table 1) . In this case the Brenner function pre- 
dicts that the formation of a dimer vacancy is energetically more prob- 
able than the formation of a single vacancy. Accordingly, for a dimer 
vacancy a considerable amount of relaxation energy is released (in par- 
ticular for the Brenner potential) during the reequilibration. 

Energy values calculated for vacancies formed on the (111) surface 
of diamond are given in Table 2. Also, schematic top views for relaxed 
vacancy configurations are depicted in Figure 2. In all cases, vacancy 
energies for the (111) surface were found to be larger than the ener- 
gies calculated for the (100) surface vacancies. In general, the Bren- 
ner function predicts somewhat smaller vacancy energies than the Ter- 
soff potential. Similar to the previous case, carbon atoms adjacent to 
the vacant site were found to be displaced outwardly upon relaxation. 

For a single vacancy created at the top layer, the displacements were 
calculated as 0.16 and 0.10 A for the three second layer atoms closest 
to the vacancy site using the Tersoff and Brenner functions, respec- 
tively. This is shown in Figure 2b. For a single vacancy formed in the 
second layer the outward displacements for the three first layer carbon 
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atoms are 0.18 and 0.09 A for the Tersoff and Brenner functions, respec- 
tively. (See Figure 2c). For a two-vacancy defect formed by removing 
two neighboring carbon atoms from the top layer, atoms adjacent to the 
vacant site were also relaxed outwardly as shown in Figure 2d. The dis- 
placements of these atoms upon relaxation are comparable to those of the 
top layer single vacancy case mentioned above. 


Conclusions 

Qualitatively speaking both potential functions produced consistent 
results which indicate that vacancies are more likely to form on the 
(100) surface than on the (111) surface. Calculations by both functions 
indicated very low formation energies for dimer vacancies on the (100) 
surface in particular. Based on energetic considerations, therefore, a 
higher concentration of vacancies is expected on the (100) surface. In 
all cases, energies calculated for single vacancies in the second layer 
were found to be higher than energies for single vacancies in the top 
layer. In general, energy differences are ^ 2 eV or higher indicating 
that a single vacancy diffusion from the top layer to the second layer is 
an energetically unfavorable process. In all cases, atoms adjacent to 
a defect site were found to relax outwards (i. e. , in the direction away 
from the vacant center) . Atoms second neighbors to a vacant site how- 
ever, exhibited only very minimal displacements during the relaxation. 

Present results may serve to build reliable models for the inter- 
pretation of experimental observations. Despite the fact that calcula- 
tions carried out here are for the low temperture limit, they are based 
on potential functions which are shown to reproduce many properties of 
diamond correctly. 
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Table 1. Calculated formation energies for vacancies on the (2x1) 
dimerized (100) surface of diamond. E Y and E v denote unrelaxed and re- 
laxed energies in eV. 




Tersoff 


Brenner 


E? 

E v 

E? 

E v 

Single Vacancy 
(Top Layer) 

2.168 

0.432 

2.686 

0.941 

Single Vacancy 
(Second Layer) 

5.062 

3.212 

4.385 

2.917 

Pair Vacancy 
(Top Layer) 

1.562 

0.483 

2.681 

0.317 
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Table 2. Calculated formation energies for vacancies on the (111) sur- 
face of diamond. E ^ and E v denote unrelaxed and relaxed energies in eV 





Tersoff Brenner 



E? 

E v 

E V v 

E v 

Single Vacancy 
(Top Layer) 

3.930 

2.939 

2.752 

2.022 

Single Vacancy 
(Second Layer) 

7.090 

5.657 

5.892 

5.189 

Pair Vacancy 
(Top Layer) 

7.924 

6.208 

5.968 

4.498 
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Figure Captions 


Figure i. Schematic top views for the (2x1) dimerized (100) surface of 
diamond, (a) defect-free surface, (6) single vacancy in the top layer, 

(c) single vacancy in the second layer and (d) dimer vacancy in the top 
layer. Large and small circles represent carbon atoms located in the top 
and second layers, respectively. 


Figure 2. Schematic top views for the (111) surface of diamond, (a) 
defect-free surface, (6) single vacancy in the top layer, (c) single va- 
cancy in the second layer and ( d ) double vacancy in the top layer. Large 
and small circles represent carbon atoms located in the top and second 
layers, respectively. Arrows indicate atomic motions in lateral direc- 
tions during relaxation. 
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Point Defect Structures and Energetics in Si Using an Empirical 
Potential 
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A recent comparative study of six different empirical potential en- 
ergy functions for silicon [1] showed that each served best for a unique 
range of application and conditions but no single potential of this set 
gave excellent fit over the broad range of applicability for which such 
potentials might be used; i.e., clusters, bulk, surfaces, defects, and 
so on. This study [1] showed that the Tersoff potentials [2,3] were ac- 
ceptable choices for the investigation of point defects and elastic 
properties, although they have very short range and do not provide per- 
fect matches with some ab initio results [1] (such as predicted struc- 
tures of Si 3 _ 6 ) . In spite of their few shortcomings, we find it is useful 
to employ such potential energy functions (PEFs) for considering atomic 
structure variations around point defects and for investigating how the 
energetics of these defects change with respect to distance from a sur- 
face or from a site of Frenkel defect formation. 

In this paper, we use the Tersoff 2 (T2) and Tersoff 3 (T3) poten- 
tials [4] to perform statics calculations (representing the low temper- 
ature limit) for the energetics of the following processes: (i) neutral 
monovacancy and divacancy formation and migration in Si; (ii) neutral 
bond- centered, site-shearing, tetrahedral, and hexagonal self- inter- 
stitial formation and migration; (iii) the variation of these energies 
with distance from a bulk site of Frenkel defect formation as well as 
from [100] and [111] surfaces. In a subsequent paper, the temperature 
dependence for these different defects will be presented where, in Monte 
Carlo simulations, one can see the structural transition from one kind 
of defect to another within a specific temperature range. Such inter- 
esting temperature-dependent behavior with its associated structural 
details is not available today from purely ab initio treatments [5-9] . 
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This is one of the definite advantages of the empirical PEF approach 
compared to present-day ab initio approaches* It is also curious to note 
that the empirical PEF technique seems to provide activation enthalpies 
for point defect formation in Si that are closer to experimental values 
than do the ab initio calculations [10] . 

Procedure 

In the present work, we have effectively confined our calculations 
to the low temperature limit by employing the Fletcher-Reeves-Polak- 
Ribiere conjugate gradient technique [11] to perform static energy min- 
imization for configurations of Si atoms as part of two kinds of static 
simulations. The first is a straightforward formation energy calcula- 
tion in which a defect is created and the atoms are then relaxed to what 
should be the minimum energy state. This relaxed energy is then compared 
to a similar system without the defect. The second is an extension to 
formation energy calculations wherein we study the process by which a 
created defect migrates from site to site. 

The formation energies per defect, A Ef, reported on here for vacan- 
cies and self-interstitials are calculated using 

A E, = E d - ^ Et (1) 

where Ed and E\> are the energies of the system with and without the de- 
fect, respectively, while and ere the total numbers of atoms in 
each system. For a vacancy, the formula assumes that the absent atom has 
moved out of the computational cell to a surface kink site where it con- 
tributes to the total binding energy without altering the energetics of 
the surface. For an interstitial, the process is reversed and an atom is 
removed from a surface kink site to an interior site within the computa- 
tional cell. 
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To calculate activation barriers for point defect migration, the 
program moves an atom along a specified path and, at each point along the 
path, it uses static minimization to calculate the relaxed energy of the 
system with the constraint that the advancing atom is (temporarily) held 
fixed. After each minimization, the other atoms return to their initial 
positions, and the moving atom advances to the next point on the path. 

Our intention is to simulate the process by which a defect is created or 
an existing defect migrates from point to point, and hence to estimate 
the potential energy barrier for this process to occur. 

The selection of the path on which to move the atom can be a study 
in itself. Ideally, the atom should experience no component of force 
perpendicular to the path before or after relaxation of the surrounding 
atoms, since a non-vanishing perpendicular force implies that there is a 
neighboring path with lower energies connecting the same two endpoints. 
This condition means that appropriate paths follow the flow lines of the 
3-D force field. Not all flow lines are appropriate, however, since the 
potential should have positive curvature perpendicular to the path so 
that it is stable. If the curvature is negative along a perpendicular 
direction, then the selected flow line effectively runs along the crest 
of a ridge, and is an unlikely candidate. Essentially, we are searching 
for a sequence of flow lines that connect a pair of stable sites through 
one or more saddle points that are energy maxima along the path and are 
stable perpendicular to the path. These saddle points are states at the 
top of the energy barriers along the path. 

The common method of doing activation barrier calculations is that 
of constrained optimization, which avoids selection of the path in ad- 
vance by moving the atom along one direction while allowing it to relax 
in the perpendicular plane. The atom then follows its own path, but not 
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necessarily one that satisfies the conditions just described. Moreover, 
the greater the relaxation of the advancing atom in the perpendicular 
plane, the more the conditions are violated. We feel it is better to 
spend more effort preselecting an appropriate path. Fortunately, the 
path required to simulate the migration of, for example, a monovacancy 
is straight and very short. Furthermore, we have found that a path con- 
structed with linear segments connecting neighboring tetrahedral sites 
is often satisfactory for longer paths in the bulk due to the short range 
of the potential and to the remaining symmetry of the crystal, even with 
a defect. On the other hand, a piece-wise linear path can be unsatis- 
factory for migration from a surface, for example, where the symmetry is 
broken. More sophisticated path selection methods are possible, but we 
have not tried them here. 

Results 

A. Point Defect Energies and Structures 

The formation energies, Ef , and relaxation displacements, S { , of 
adjacent atoms for several point defects axe given in Table I. Fig- 
ures 1 and 2 provide representations of some different local atomic 
configurations encountered. The tetrahedral (T) and split (S) mono- 
vacancy structures plus the stable (D) and metastable (DO divacancy 
structures are given in Figure 1. During relaxation, the atoms neigh- 
boring the T-vacancy move outwaxds while they move inwards for the 
S-vacancy. We attribute this difference to the presence of a cen- 
tral atom in the split form that pulls surrounding atoms inward to 
help fill the void. We see from Table I that the S-monovacancy has 
the smaller total energy at the expense of increased atomic dis- 
placement . For the different interstitial species, Figure 2 illus- 



trates the different types of displacements, Si , listed in Table I. 


B. Point Defect Migration 

Figure 3a illustrates the path used to simulate the migration of a 
monovacancy . Figures 4a and 4b represent the simulation results for the 
T2 and T3 PEFs, respectively. From 0 to 1 on the abscissa corresponds to 
a displacement of 2.35 A (from one tetrahedral site to another) . Note 
that the vacancy structure switches discontinuously from the tetra- 
hedral structure to a split form near mid step. For the T-vacancy, the 
moving atom initially has three neighbors, while for the S-vacancy it 
has six. It is precisely when the moving atom encounters the additional 
neighbors that a jump in energy and a change in structure occurs. 

In this case, such behavior is a probable consequence of the short 
range and sharp cutoff function characteristics of the Tersoff poten- 
tials, as well as the limitations of the static minimization procedure. 
For example, if the relaxation of a configuration of atoms to its low- 
est energy state requires more than a small motion of any of the atoms, 
then the conjugate gradient procedure (as well as any other multidimen- 
sional static minimization technique) is likely to relax the atoms into 
a metastable state, since it cannot "see" the lower energy state. Other 
problems arise from the form of the cutoff function in the transition 
region of width fi from R c — \i to R c Although the function is smooth, its 
transition to zero is very sharp and its slope has corners at the end- 
points of this region (i.e., the second derivative is discontinuous), 
which can be problematic for routines using the gradient directly or in- 
directly by quadratic interpolation. The cutoff function may be modi- 
fied to eliminate this behavior, but we have not done this for the work 
presented here, tfe have instead monitored the number of neighbors of a 
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defect atom (for example, the moving atom in a migration calculation) 
both before and after relaxation of the lattice. If the number of neigh- 
bors changes, either during relaxation or as a migrating atom moves, 
then there is the possibility that the energy and structure will exhibit 
large changes due to the nature of the cutoff and the limitations of the 
simulation method. 

In spite of the difficulty in providing an unambiguous interpreta- 
tion, the migration barriers are found to be 0.25 eV and 1.60 eV for T- 
and S-vacancies, respectively, with the T2 PEF, and 2.5 eV and 2.7 eV for 
T- and S- vacancies, respectively, with the T3 PEF. (These values were 
obtained from the energy versus position data generated by the program.) 
The energies along the migration path are higher for the T3 potential 
than for the T2 potential because T2 underestimates bond-bending forces 
(the shear modulus computed for a single crystal with the T2 potential 
is an order of magnitude lower than the accepted value) , while T3 gives a 
much better match [2,3,12,13] . 

Figure 3b shows the process used to model the divacancy migration. 

An atom adjacent to one of the vacant T-sites moves straight to the near- 
est open site, which forms a metastable divacancy, and then continues to 
the next vacant site for a net shift of the divacancy center. This de- 
scribed path may not be a true flow line with vanishing perpendicular 
force; however, only small changes in energy were observed when attempt- 
ing to refine the path by using constrained optimization. The migration 
energies for divacancies are shown in Figures 5a said 5b for the T2 PEF 
said the T3 PEF, respectively. Note that the migration is a two-step pro- 
cess and that each step is similar to the migration of a monovacancy. 

The activation barriers for divacancy migration were found to be 2.2 eV 
for the T2 PEF and 3.6 eV for the T3 PEF using T-vacancies. 
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The migration of a self-interstitial atom from one tetrahedral site 
to another was simulated yielding the results given in Figure 6a for 
the T2 PEF and Figure 6b for the T3 PEF , respectively. The calculation 
includes all three of the bond-line interstitial types since the nom- 
inal H*- and H-sites are on the line connecting adjacent tetrahedral 
sites. Although milder, the discontinuities in the relaxed energy axe 
due to the reasons already discussed for the vacancy defect. Gener- 
ally, we cannot unambiguously favor one portion of the curve over an- 
other, but in the case of a discontinuity as narrow as the pedestal in 
Figure 6a (less than 0.02 A) not only would a classical simulation at 
temperatures greater than zero wash it out completely, but also the 
atom will not be confined to such a small region quantum mechanically 
since the corresponding uncertainty in the momentum of the atom will be 
large. If we eliminate the sharp pedestal, the revised estimate of the 
T-interstitial formation energy is 4.0 eV for T2, although this struc- 
ture is no longer expected to be stable. The migration barriers are 0.5 
eV for the H*-interstitial for T2, and 1 . 1 eV for the T-interstitial for 
T3. Other barriers are listed in Table II. 

Figure 7 shows how the program simulates the migration for the S- 
interstitial. The atom labeled n l" moves into the vacant T-site (shown 
as a thin, white circle), which causes the other atom in the split in- 
terstitial (labeled M 2 M ) to move to a nearby T- or H*-interstitial po- 
sition. The S- interstitial is thus converted to one of the other inter- 
stitials. This is one half of the process; the other half of the process 
is the analogous inverse process that converts the simple interstitial 
into an S-interstitial at an adjacent location. The activation barriers 
for the total process are given in Table II. 


The energy profiles for point defect movement inwards from the sur- 



fac© for the two PEFs and the two selected surface orientations are 
given in Figure 8 for T-vacancies and Figure 9 for T-interstitials . We 
see that these profiles are fully stabilized at roughly 3 atomic dis- 
tances from the surface. There is evidence of low energy states in the 
surface layers where defects would, therefore, tend to concentrate. 

The [100] , T2 calculations actually exhibit negative vacancy and in- 
terstitial formation energies (clipped on the graphs) near the surface 
layer, which implies that the surface is ripe for reconstruction using 
these parameters . Finally, the energy profile for the T- vacancy and T— 
interstitial formation as a function of distance from a Frenkel defect 
formation site is presented in Figure 10. The activation energy for for- 
mation of this type of defect is 1.0 eV for T3 but only 0.38 eV for T2. 

The intricate detail in the curves in Figures 8-10 is the result of 
the short range of the PEFs. As an atom moves along the bond line from 
a tetrahedral site, it first comes under direct influence of the atoms 
in the hexagonal ring surrounding the H-site, then it flips through the 
ring (i.e. , it moves across the H-site), and continues to the next T- 
site. Since the potentials that we used have short ranges, the number 
of neighboring atoms (as well as which atoms are the neighboring atoms) 
that contribute to the energy changes signif icantly ; a potential with 
greater range would wash out some of this detail by averaging into the 
calculation significant contributions from neighbors beyond the first 
shell. 

In the development of both PEFs, cohesive energy data for real and 
hypothetical bulk structures, the bulk modulus, and the bond length 
in the diamond structure were used [2,3] . However, the T3 results are 
preferable to the T2 results, because the T3 PEF satisfies the added 
constraint that it reproduce all three elastic constants of bulk sili- 


con to within 20'/,. This greatly improves the shear modulus, which is too 
low by an order of magnitude using T2 . 

For an example calculation with the T3 results. Table I shows that 
the S-vacancy and the T-interstitial have the lowest formation energies 
for T — > OK. Assuming that the formation entropy for all point defects is 
A Sf ~ 5k B , we have 


C*(T ) = N a exp(ASf/kB)exp(-AHf/k B T) (2a) 

Cy 0 (T) = 7.5 x 10 2i exp(—3.52eV/kBT)cm~ 3 =$■ Cy 0 (TM ) = 2.78 x 10 14 cm -3 (26) 

Cjo[T) = 7.5 x 10 24 erp(— 3.48eVy&BT)cm -3 => C}o(T M ) = 3.75 x 10 14 cm- 3 (2c) 

for the equilibrium vacancy and interstitial concentrations, respec- 
tively, as functions of temperature and at the melting point, T m = 1683, 
and where N a = 5 x 10 22 cm ~ 3 . 

The actual picture is a little different from this because charged 
vacancy species must also be taken into account. The energy levels in 
the silicon band gap for these species are known [14] and it can be shown 
that at the Si melting temperature, the total vacancy content, Cy^ tot ), 
is given by [15] 

Cv(tot)( T M) = Cyo + C'y- + C*y= + C'y+ w 20 C'yo{T M ) (3) 


If there are analogous levels in the band gap for the charged inter- 
stitial states, these are not known at the present, so we shall assume 
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that equations 2b represent the total interstitial content at Tm • In 
such a case, Frenkel defect annihilation for T < Tm will occur be- 
tween ~ 5 x 10 15 cm -3 total vacancies and ~ 3.75 x 10 14 cm -3 neutral 
interstitials grown-in at T — Tm leaving a net vacancy excess of 
~ 4.6 x 10 15 cm ~ 3 vacancies at lower temperatures. Of course, calcula- 
tions of the thermochemical properties as a function of temperature will 
yield slightly different results, but these numbers might be considered 
to be "ball park" estimates. 

With respect to defect migration using T3, the T-vacancy is the 
most stable monovacancy with a migration energy of 2.5 eV, while the S- 
interstitial variety is a faster moving species with an activation en- 
ergy of 0.27 eV, although its equilibrium population will be appreciably 
lower than the T-interstitial species. We must await higher temperature 
calculations to determine which of these defect types maintain entropi- 
cally distinguishable features, and thus which are stable defects that 
are present in the lattice and actively involved in transport . 

IV. Conclusions 

It is important to recognize that these 0 K energy barriers, as cal- 
culated by the present procedures, cannot be expected to correspond al- 
ways to what high temperature simulations calculate as the activation 
energy for the process. It is this latter quantity that will be acces- 
sible experimentally. The prime value of the present work thus lies in 
its internal comparison between different defects rather than its match 
with experimental data. 

Preferring the T3 PEF for point defect calculations to the T2 PEF, 
the S- and T-monovacancies are found to have comparable A Ef + A E m 
values with A Ej v = 3.70 eV, A E^ v — 3.52 eV, AE^ V = 2.50 eV, and 
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A v = 2.70 eV . The main distinguishing feature between these two 
species is that the first neighbor shell of atoms relaxes outwards for 
the T-vacancy and inwards for the S-vacancy. A somewhat similar sit- 
uation exists for the T-, H*-, and S- interstitials with A Ejj = 3.48 
eV, A Efj = 4.08 eV, A E s f I = 4.41 eV, A j = 1.1 eV, AE%* T = 0.56 
eV ernd Aj r = 0.27 eV. We find the H-interstitial to be unstable 

0 eV) and the first neighbor shell of atoms relaxes outwards 
for all species. It was also found that the activation energy for T-T 
Frenkel defect formation is 1.0 eV and that the variation of the above 
energies occurs only within ~2-3 atomic distances of a [111] or [100] 
surface. 
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Table I: This table shows the defect formation energies, A Ef f defined 
by equation 1, and the relaxation displacements of the atoms near the 
defect, Si, (in order i=l,2, . . . corresponding to those illustrated in 
Figures 1 and 2) . The displacements are given relative to the cubic lat- 
tice parameter, a = 5.43 A. A negative displacement indicates relaxation 
inward (i.e., toward the defect), We use T, H, H*, and S to refer to 
tetrahedral, hexagonal, bond-centered, and split forms, respectively. 

D and D J refer to the stable and metastable divacancies. We designate 
the Frenkel defects according to "vacancy type + interstitial type". 

Table II: Activation energy data for migration of point defects. The en- 
ergy differences were obtained from the numerical values calculated by 
the migration simulation program. We have given two energies and a se- 
quence of sites describing the corresponding migration paths for the H*- 
interstitial and S-interstitial defects. Some of the defects are unsta- 
ble with respect to migration, which suggests that these configurations 
are multidimensional saddle points rather than true energy minima. 
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Defect 

Tersoflf 2 
A £>(eV) S/a 

A£ r (eV) 

Tersoflf 3 
<5/a 

vacancies 





T 

2.81 

0.0043 

3.70 

0.045 

S 

1.46 

-0.080 

3.52 > 

-0.050 

divacancies 





D 

4.22 

0.0043 

5.57 

0.044 

D' 

6.22 

<5,/a=0.0049 

7.64 

<5,/o=0.045 



6/o= 0.0049 


<5 2 /o=0.03 1 

interstitials 





T 

5.04 (4.0) 

6/a= 0.032 

3.48 

6/a=(0) 



6 2 /a=(0) 


<5 2 /o=0.045 

H 

3.68 

0.032 

4.63 

0.043 

H* 

3.50 

6/0= 0.036 

4.08 

<5,/tf=0.040 



6 2 ta= 0.021 


6/ a=Q .021 



<5^0=0.013 


6/a=0.0\9 

S 

2.46 

<5,/a=0.082 

4.41 

( 5/a=0.063 



S 7 la= 0.074 


6/ a=Q. 04$ 



<5-,/o=0.019 


6/a=0.021 

Frenkel 





T+T 

7.85 


7.19 


T-H 

6.49 


8.34 


T+H* 

6.31 


7.78 
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Table H 


Point Defect 
Type 

Tersoff2 
A E m (eV) 

Tersoff 3 
A E m (eV) 

vacancies 

T 

0.25 

2.5 

S 

1.60 

2.7 

divacancies 

D 

2.2 

3.6 

D' 

0.21 

1.55 

interstitials 

T 

unstable 

1.1 

H 

unstable 

unstable 

H* 

0.5 (H*->T->H*) 

0.56 (H*-»H-»H*) 


0.17 (H*— »H— »H*) 

0.02 

S 

1.5 (S-4H*^T->H*->S) 
1.2 (S— »H*— »S) 

0.27 (S-VT-*S) 



Figure Captions 
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Figure 1: This figure shows the four different vacancy configurations: 
(a) tetrahedral and (b) split monovacancy structures, and the (c) sta- 
ble and (d) metastable divacancy structures, designated T, S, D, and D’, 
respectively. The atoms around the tetrahedral vacancy relax outwards 
from their initial locations by a distance 8 = 0.0043a (a = 5.43 A is the 
cubic lattice parameter) for T2 and by 8 = 0.045a for T3. In the split 
structure, one of the atoms next to the vacant tetrahedral site moves 
halfway towards the empty site, thus splitting the vacancy between two 
sites. The surrounding atoms relax inwards toward the site opposite the 
intermediate atom by 8 = 0.080a for T2 and by 8 = 0.050a for T3. The di- 
vacancy formation energies are largely consistent with the number of 
bonds broken, and the relaxation of the atoms surrounding the divacancy 
is similar to that observed for T-vacancies. The atoms around the D’- 
divacancy relax by different amounts, so the figure shows two displace- 
ments, and 82 . 

Figure 2: Four interstitial structures cure shown. Figure (a) shows the 
tetrahedral (T) structure, (b) shows the hexagonal (H) structure, (c) 
shows the alternate bond-centered or axial hexagonal (H*) structure, 
and (d) shows the split-site form. The interstitial atoms in (a)-(c) 
are white and marked with an "I", while in (d) the two atoms that form 
the split structure are white and marked with "S" . The figures also show 
how sets of atoms around the defect relax. Each set consists of atoms 
that relax with the same size of displacement (shown as vectors with 
length 8 i ) . When we show more them one set, we have ranked them in de- 
creasing order of displacement magnitude emd marked the members of each 
set with a number corresponding to their rank. The H*- and S- intersti- 
tial structures have relaxation displacements that are significant for 
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more atoms than for the T- and H- interstitials, so we have indicated up 
to three different displacements* Note that for the T2 potential, only 
the nearest neighbors undergo any appreciable relaxation around a T- 
interstitial, while for the T3 potential, these atoms move very little 
and the atoms in the second neighbor shell move outwards. 

Figure 3: Paths selected to simulate (a) vacancy and (b) divacancy mi- 
gration. 

Figure 4: The two graphs show the results of the vacancy migration cal- 
culations for bulk silicon. We present the vacancy energy both before 
and after relaxation as it migrates from one T-site to an adjacent one 
for (a) T2 and (b) T3. Note that the vacancy switches to the split form 
in mid step. 

Figure 5: The two graphs show the energy as a divacancy migrates in bulk 
silicon via the process illustrated in Figure 3b. The energy before and 
after relaxation is shown for (a) T2 and (b) T3. 

Figure 6: The graphs show the interstitial energy as the defect migrates 
from one tetrahedral site to an adjacent one in bulk silicon. Note that 
the H*- and H-sites lie along this path. The unrelaxed and relaxed ener- 
gies are presented for (a) T2 and (b) T3. 

Figure 7: The arrows illustrate the first half of the S-interstitial 
migration process: the atom labeled "i" moves into the vacant T-site 
(shown as a white circle), which causes the other atom comprising the 
split interstitial (labeled M 2 n ) to move to a nearby interstitial po- 
sition (such as T or H*) , converting the S-interstitial to one of the 
other types. The second half of the migration step is an analogous in- 
verse process that converts the simple interstitial into a split inter- 
stitial about a nearby T-site. 


49 



/ 


Figure 8: The graphs show the vacancy migration energies as the defect 
migrates inwards from [111] and [100] type surfaces both before relax- 
ation (ur) and after relaxation (r) of the surrounding atoms. Illus- 
trated are migration from [111] for T2(a) and T3(b), and migration from 
[100] for T2(c) and T3(d). The steps indicated on the horizontal axes 
correspond to completed hops from tetrahedral site to tetrahedral site. 
Each site is successively deeper into material. 

Figure 9: The graphs show the interstitial migration energies as the 
defect migrates inwards from [111] and [100] type surfaces both before 
relaxation (ur) and after relaxation (r) of the surrounding atoms. Il- 
lustrated are migration from [111] for T2 (a) and T3 (b) ; and migration 
from [100] for T2 (c) and T3 (d) . The steps indicated on the horizontal 
axes correspond to completed hops from tetrahedral site to tetrahedral 
site. Each site is successively deeper into the material. 

Figure 10: The graphs illustrate the formation of the T-vacancy/T- 
interstitial Frenkel defect for (a) T2 and (b) T3. tfe show both relaxed 
and unrelaxed curves . 
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